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Introduction 


This  report  describes  the  progress  and  accomplishments  for  the  period  June  1,  2002  to 
August  31,  2004  for  the  MURI  grant  “Accurate  Theoretical  Predictions  of  the  Properties  of 
Energetic  Materials”  (No.  DAAD 19-02- 1-0 176).  This  is  a  multi-university,  comprehensive 
theoretical/computational  research  program  to  develop,  validate,  benchmark,  and  apply  methods 
and  models  that  will  provide  predictive  capabilities  for  energetic  materials.  The  thrust  of  the 
work  is  the  development  of  atomic-level  models  and  ab  initio  quantum  chemistry  methods  that 
are  generally  applicable  to  the  chemical  decomposition  of  condensed-phase  energetic  materials 
under  extreme  conditions. 

The  approaches  include  quantum  mechanics,  molecular  modeling,  Monte  Carlo,  and 
molecular  dynamics,  to  yield  state-of-the-art  methods  specifically  designed  for  and  tailored  to 
target  DoD  energetic  materials  research  needs.  The  following  is  a  list  of  the  principal 
investigators,  universities,  and  topics  that  make  up  the  overall  MURI  project: 

•  University  of  Missouri-Columbia;  Professor  John  E.  Adams;  Gas-liquid  heat  transfer, 

•  University  of  Maryland;  Professor  Herman  E.  Ammon;  Crystal  models 

•  University  of  Florida;  Professor  R.  Bartlett;  Quantum  Chemistry 

•  North  Carolina  State  University;  Professor  Donald  W.  Brenner;  Reaction  potentials 

•  University  of  Illinois;  Professors  David  Ceperley  and  Richard  Martin;  Solid-state  quantum 
mechanics 

•  University  of  Missouri-Columbia,  Professor  Donald  E.  Thompson  (PI);  Reactions  dynamics 
and  simulations 

•  University  of  Minnesota;  Professors  Donald  G.  Truhlar  and  Christopher  J.  Cramer;  Solvation 
models 

The  focus  of  the  project  is  on  developing  accurate  methods  for  simulating  physical  and 
chemical  processes  in  condensed-phase  energetic  materials.  The  MURI  team  is  working  closely 
with  researchers  at  DoD  labs,  which  are  contributing  to  the  theoretical  efforts,  providing  data  for 
testing  of  the  models,  or  aiding  in  the  transition  of  the  methods,  models,  and  results  to  DoD 
applications.  The  major  goals  of  the  project  are: 

•  Models  that  describe  phase  transitions  and  chemical  reactions  in  condensed  materials. 

•  Ab  initio  predictions  of  structures  and  properties  of  solids  at  high  temperatures  and  pressures. 

•  Methods  to  predict  mechanical  properties  and  physical  changes  in  condensed  phases,  including 
thermodynamic  quantities,  phase  diagrams,  vibrational  frequencies,  tensile  strength,  stress/strain 
relations,  and  densities. 

•  Simulation  methods  to  predict  chemical  decomposition  in  condensed  phases,  particularly 
ignition  and  sensitivity  in  response  to  heating  and  shocking. 

•  Methods  for  predicting  temperatures  of  the  condensed  phases  and  flames  resulting  from 
physical  and  chemical  changes,  including  a  predictive  model  for  the  “heat  feedback”  from  the 
flame  to  burning  surface. 

•  Methods  for  predicting  solvation  and  separation  for  energetic  materials  in  supercritical  fluids. 

Technical  Progress 

The  following  are  descriptions  of  the  components  that  make  up  the  overall  project. 
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Structure  and  Density  Predictions  for  Energetic  Materials 

(Herman  L.  Ammon,  Department  of  Chemistry  and  Bioehemistry,  University  of  Maryland, 
College  Park,  MD) 

The  diseovery  of  new  energetie  materials  ean  be  faeilitated  with  eomputer  modeling  and 
simulation  teehnology  for  the  identifieation  of  eompounds  with  signifieant  advantages  over 
materials  eurrently  in  use.  The  quantitative  estimation  of  properties  sueh  as  the  heat  of 
formation,  density,  detonation  veloeity,  detonation  pressure  and  sensitivity  ean  sereen  potential 
energetie  candidates  and  permit  the  selection  of  only  the  most  promising  substances  for 
laboratory  synthesis,  measurement  of  properties,  scale-up,  testing,  etc.  Thus  far,  our  work 
primarily  has  involved  developing  computer  codes  and  procedures  for  the  prediction  of  the 
crystal  structures  of  small,  covalent  C-H-N-O-F  containing  molecules.  The  calculations  give 
solid-state  densities  and  heats  of  formation  and  the  hypothetical  structures  can  provide  valuable 
information  of  sensitivity  estimates. 

Our  model-MOLPAK-refinement  procedure  to  predict  the  most  probable  crystal  structure 
from  a  rigid  molecular  model  involves  three  steps:  (1)  a  model  for  the  compound  of  interest  (the 
search  probe)  is  obtained  from  ab  initio  quantum  mechanical  geometry  optimization  (usually 
B3LYP/631G*  or  63 IG**);  (2)  determination  of  thousands  of  possible  crystal  structures  for  the 
search  probe  (MOLPAK  program,  MOLecular  PAcKing);  (3)  refinement  (WMIN,  monopole 
electrostatics,  or  DMAREL,  distributed  multipole  electrostatics)  of  the  unit  cell  parameters, 
search  probe  orientation  and  position  by  lattice  energy  minimization  for  the  best  of  the  crystal 
structures  derived  in  step  2.  The  most  probable  structure  is  identified  by  a  combination  of 
(lowest)  crystal  lattice  energy  and  (highest)  density.  In  tests  of  molecules  with  known  crystal 
structures,  the  experimental  structure  is  identified  by  the  energy  criterion  in  80%  of  the  cases 
and,  for  the  other  20%,  the  correct  structure  is  among  the  half-dozen  lowest  energy  predicted 
structures. 

In  WMIN,  the  crystal  lattice  energy  (E)  is  calculated  from  the  sum  of  all  atom-to-atom 
interactions  between  a  central  molecule  and  its  neighbors. 

Ewmin  =  X{332.17[qiqjrij'^]  -  AAiPj'^  +  BiBjexp[-(Ci  +  Cj)rij]} 

The  terms  are:  qi  =  Coulombic  charge  on  atom  i;  pj  =  atom  i  to  j  distance;  Aij  =  (Ai*Aj)  ,  Bij  = 
(Bi*Bj)^^^,  Ai  and  Bi  are  empirical  coefficients  for  atom  i;  Ci  is  similar  to  a  van  der  Waals  radius 
for  atom  i.  The  DMAREL  potential  is  similar  with  a  distributed  multipole  derived  electrostatic 
term  (63 IG**  basis  set)  in  place  of  the  WMIN  monopole.  The  A  and  B  coefficients  are 
optimized  by  fitting  experimental  and  calculated  crystal  lattice  data  with  simplex,  gradient,  and 
least  squares  methods.  Despite  the  paucity  of  experimental  enthalpies  of  sublimation  to  provide 
lattice  energy  standards,  the  parameterizations  provide  a  good  set  of  coefficients,  which  translate 
into  more  sensitive  and  definitive  lattice  energies,  structure  refinements  and  accurate  predictions. 
Currently,  A  and  B  coefficients  have  been  determined  for  69  C,  H,  N,  O,  F,  S  and  Br  atom  types. 

The  WMIN  program,  developed  at  Oak  Ridge  National  Laboratory  about  1980,  is 
undergoing  a  complete  re-write.  The  new  code  will  incorporate  modular,  modem  Fortran  and 
have  improved  transparency  and  portability.  Major  additions  will  include  A  and  B  cross  terms 
for  all  atom-to-atom  interactions  to  facilitate  the  handling  of  special  contacts  such  as  H-bonding, 
analytical  derivatives  in  place  of  the  present  numerical  quantities,  facilitated  handling  of  multiple 
non-covalently  linked  groups  in  the  unit  cell  and  intramolecular  optimization. 
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The  search  probe  model  obtained  from  ab  initio  MO  geometry  optimization  describes  a 
molecule  free  of  crystal  packing  effects.  Although  a  rigid  probe  is  adequate  for  many  packing 
calculations,  there  are  cases  in  which  the  incorporation  of  conformational  flexibility  as  part  of 
the  packing  process  is  required.  The  flexibility  problem  has  been  addressed  by  the  development 
of  the  ROTPAK  program  (ROTational  PAcKing)  that  allows  defined  conformational  changes  to 
accompany  the  packing  processes.  The  “best”  structures  represent  minima  in  the  total  (intra  plus 
intermolecular)  energy.  Although  conceptually  simple,  the  process  is  complicated  by  (1)  the 
determination  of  the  intramolecular  energy  and  (2)  achieving  a  proper  balance  between  the  intra 
and  intermolecular  energies.  With  the  use  of  the  PM3  semi-empirical  method  to  estimate 
intramolecular  energies,  ROTPAK  has  been  used  to  successfully  predict  the  crystal  structures  of 
the  nitrocubanes,  l,3,5-triazido-2,4,6-trmitrobenzene,  N-cyano-N-nitrotolylamine  and  2-methyl- 
4,5-dmitroacetamide  starting  with  B3LYP/631G*  optimized  models. 

Sensitivity,  the  ease  with  which  a  material  undergoes  a  violent  reaction  or  detonation,  can 
be  triggered  by  numerous  stimuli  such  as  impact,  shock,  friction,  thermal  and  electrostatic 
sources  of  energy.  A  number  of  molecular  structure-sensitivity  relationships  have  been  found 
over  the  years  and  good  correlations  observed  within  families  of  energetic  materials.  Sensitivity 
is  perhaps  the  most  complicated  and  least  well  understood  of  the  various  derivative  properties  of 
energetic  materials.  The  following  areas  will  be  investigated.  (1)  We  will  further  explore  the 
utility  of  the  density  of  states  relationship  of  Kunz,  which  requires  the  crystal  structure  and 
appropriate  ab  initio  crystal  lattice  calculations.  (2)  Additionally,  possible  relationships  between 
impact  and  friction  sensitivity  and  bond  strength  plus  lattice  energy  will  be  investigated.  The 
basic  idea  is  that  impact  has  an  important  frictional  or  heat  component  that  initially  causes 
disruption  (melting  or  deformation)  of  the  crystal  lattice  followed  by  homolytic  cleavage  of  the 
weakest  bond  in  the  molecule  followed  by  detonation.  Lattice  energies  are  available  from 
WMIN/DMAREL  calculations  and  bond  energies  (e.g.  C-NO2,  N-NO2,  O-NO2)  from  ab  initio 
calculations  of  the  energies  of  the  appropriate  framework  and  nitro  group  radicals.  (3)  We  have 
observed  previously  a  correlation  between  the  impact  sensitivity  and  hydrostatic  compressibility 
(from  molecular  dynamics  calculations)  for  a  variety  of  energetics.  For  example,  the  impact 
sensitive  PETN  is  more  compressible  than  the  relatively  insensitive  TATB.  This  will  be 
investigated  further.  (4)  The  lattice  potentials  developed  can  be  utilized  to  explore  the 
relationship  between  crystal  orientation  and  detonation.  In  PETN,  for  example,  the  orientation  to 
shock  initiation  and  detonation  is  consistent  with  steric  hindrance  to  sheer  at  the 
molecular/crystal  level. 

Continuing  and  future  work 

1.  MOLPAK  and  ROTPAK  preliminary  searches.  Adjust  the  intermolecular  coefficients  to 
provide  better  hypothetical  structures  for  subsequent  lattice  optimization. 

2.  ROTPAK.  Continue  development  with  near-term  focus  on  intramolecular  energy  evaluation 
and  handling  of  multiple-bond  flexibility. 

3.  MOLPAK/ROTPAK.  Merge  rigid  and  flexible  molecule  codes  to  a  single  unit  with  complete 
capabilities  of  each. 

4.  Continue  force  field  coefficient  optimization  for  WMIN  and  DMAREE.  Extend  to  new 
functional  groups  and  H-bonding. 

5 .  Structure  selection.  Investigate  other  criteria  (presently  use  lattice  energy  and  density),  such 
as  a  comparison  of  patterns  of  intermolecular  contacts  with  known  crystal  structures,  to  identify 
the  best  (correct)  hypothetical  structure  from  a  prediction  set. 
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6.  Complete  the  reeoding  of  WMIN  with  a  new  lattiee  refinement  program  with  the  following 
properties:  (1)  Fortran-90  eode  that  will  exeeute  on  any  eomputer  platform;  (2)  analytieal 
derivatives;  (3)  eoupled  inter  and  intramoleeular  refinement;  (4)  separate  atom(i)-atom(j) 
intermoleeular  potentials;  and  (5)  anisotropie  potentials. 

7.  Predietion  paekage.  Couple  the  various  proeedures  into  a  seamless  paekage  that  would 
inelude:  (1)  strueture  optimization;  (2)  strueture  predietion  and  seleetion;  (3)  solid-state  heat  of 
formation  oaleulation;  and  (4)  property  ealeulations  (e.g.  Cheetah). 

8.  Sensitivity.  The  strueture  predietion  and  lattiee  potential  work  will  serve  as  a  platform  to 
examine  impaet/shoek  and  frietion  sensitivity.  Several  meehanisms  that  will  be  investigated  are 
eompressibility,  (weakest)  bond  breaking  and  lattiee  energy  and  sterie  hindranee  to  sheer. 

Personnel: 

Herman  Ammon,  Professor 
Zuyue  Du,  researeh  assoeiate 

Sayta  Prasad,  sabbatieal  assoeiate  (Physies,  Ranehi  Univ,  India) 

Ed  Wells,  graduate  student  assoeiate  and  programmer 
Nieole  Dueker,  undergraduate  assistant 
Interactions  with  the  energetic  material  community: 

Dr.  Betsy  Riee  (Aberdeen  Researeh  Laboratory) 

Dr.  Ed  Byrd  (Aberdeen  Researeh  Laboratory) 

Dr.  Alfred  Stern  (Senior  Researeh  Seientist/Teehnieal  Consultant  for  Energetie  Materials, 
NAVSEA,  Indian  Head,  MD) 

Dr.  Wiliam  Koppes  (NAVSEA,  Indian  Head,  MD) 

Dr.  Rao  Surapaneni  (Chief,  Explosives  Researeh  and  Teehnology,  Pieatinny  Arsenal) 

Dr.  Paritosh  Dave  (Pieatinny  Arsenal,  NJ) 

Dr.  Robert  Chapman  (Naval  Air  Warfare  Center,  China  Lake,  CA) 

Dr.  Philip  Eaton  (Univ.  of  Chieago) 

Dr.  Harold  Sheehter  (The  Ohio  State  Univ.) 

Puhlications: 

•  “Strueture  and  Density  Predietions  for  Energetie  Materials,”  J.  R.  Holden,  Z.  Du  and  H.  L. 
Ammon,  in  Energetic  Materials.  Part  1:  Decomposition,  Crystal  and  Molecular  Properties,  P. 
Politzer  and  J.  S.  Murray,  eds.,  pp.  183-213,  Elsevier,  Amsterdam,  2003. 

•  “Crystal  Strueture  Predietion  of  Small  Organie  Moleeules:  a  Third  Blind  Test,”  G.  Day,  W  .D. 
S.  Motherwell,  H.  L.  Ammon,  J.  D.  Dunitz,  A.  Dzyabehenko,  P.  Erk,  A.  Gavezzotti,  D.  W.  M. 
Hofmann,  E.  J.  J.  Leusen,  J.  P.M.  Lommerse,  W.  T.M.  Mooij,  S.  L.  Priee,  H.  Seheraga,  B. 
Sehweizer,  M.  U.  Sehmidt,  B.  P.  van  Eijek,  P.  Verwer,  Acta  Crystallogr.,  2004,  in  preparation. 


First  Principles  Simulations  of  Energetic  Materials 

(David  Ceperley  and  Richard  Martin,  Department  of  Physics,  University  of  Illinois,  Urbana- 
Champaign) 

The  overall  goal  of  the  group  at  the  University  of  Illinois  Urbana-Champaign  is  to 
enable  first-principles  simulations  of  energetic  materials  that  are  accurate  and  useful.  These 
methods  will  become  robust  tools  for  prediction  of  the  properties  of  materials  under  extreme 
conditions,  especially  high  temperature  and  pressure.  The  work  focuses  primarily  upon 
quantum  Monte  Carlo  (QMC)  methods,  the  only  known  approaches  that  can  simulate  materials 
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directly  from  the  fundamental  equations  for  interacting  nuclei  and  electrons.  Two  QMC 
methods  are  under  development:  Path  Integral  Monte  Carlo  and  Coupled  Electron  Ion  Monte 
Carlo.  The  primary  applications  are  to  materials  containing  first  row  elements  H,  C,  N,  and  O, 
including  selected  benchmark  calculations  for  definitive  comparison  with  experiments,  and 
simulations  of  materials  under  conditions  beyond  the  capabilities  of  direct  experimental 
analysis. 

Code  &  Algorithm  Development: 

A  particular  focus  in  our  work  under  this  grant  is  to  develop  new  methods  and  codes. 
During  this  year  we  have  made  significant  progress  in  writing  new  quantum  Monte  Carlo  codes 
that  will  allow  calculations  on  much  more  complicated  systems  and  to  effectively  use  more 
powerful  computers.  Part  of  this  work  involves  the  PIMC  codes  we  have  developed  to  simulate 
systems  of  electrons  and  protons  at  high  pressures  and  temperatures.  This  is  a  uniquely  powerful 
tool  for  simulations  as  demonstrated  by  our  work  making  detailed  comparisons[l]  with  other 
types  of  models  and  predictions  of  the  results  of  shock  wave  experiments  on  deuterium. 

There  is  a  need  to  develop  new  approaches  that  allow  the  simulations  to  be  feasible  at 
lower  temperatures  and  for  materials  other  than  hydrogen.  Our  recent  progress  is  in  the 
development  of  the  coupled  electronic-ionic  Monte  Carlo  (CEIMC)  simulations. [3]  In  this 
method,  one  moves  the  ions  classically  while  calculating  the  electronic  potential  (i.e.  Born- 
Oppenheimer)  with  QMC.  In  the  last  2  years,  this  method  has  been  tested  on  molecular  and 
metallic  hydrogen  composed  of  up  to  100  atoms  at  temperatures  as  low  as  300K.  Our  studies  of 
the  method  applied  to  hydrogen,  show  that  it  has  the  same  order  of  computational  demands,  but 
can  be  more  accurate  than  Car-Parrinello  plane-wave  methods.  The  processing  power  of  current 
multi-processors  is  enough  that  significant  applications  can  already  be  envisioned.  We  have 
developed  “reptation”  algorithms  allowing  us  to  determine  energy  differences  and  forces  much 
more  accurately  than  previously  and  new  more  efficient  methods  of  moving  the  electrons.  We 
have  also  implemented  (for  the  first  time)  constant  pressure  QMC  methods,  much  more 
convenient  for  simulation  of  energetic  materials.  Tests  for  non-hydrogenic  systems  are  needed  to 
find  the  performance  of  the  algorithms  on  a  broader  spectrum  of  applications.  We  have  begun 
constructions  of  such  methods  to  be  able  to  apply  the  CIEMC  method  to  energetic  materials  at 
finite  temperatures,  in  particular  to  allow  use  of  more  general  wavefunctions  and 
pseudopotentials. 

The  calculation  of  forces  is  an  outstanding  problem  in  Monte  Carlo  methods.  One  of  the 
main  reasons  keeping  QMC  methods  from  becoming  a  realistic  alternative  to  density  functional 
theory  (DET)  based  techniques  is  the  lack  of  an  efficient  estimator  for  nuclear  forces:  e.g.  they 
are  needed  in  molecular  dynamics  simulations.  In  the  framework  of  independent-particle  DET 
theories,  nuclear  forces  are  computed  with  the  Hellman-Eeynman  theorem:  the  force  acting  on  a 
given  nucleus  is  equal  to  the  electrostatic  force  due  to  the  charge  distribution  surrounding  it.  This 
simple  result  cannot  be  systematically  exploited  in  QMC  due  to  the  divergent  variance  of  the 
corresponding  estimator.  Possible  solutions  to  the  QMC  force  problem  have  been  suggested  but 
none  of  them  has  gone  beyond  applications  to  diatomic  systems.  We  have  come  up  with  two  new 
estimators  [4]  based  on  the  observation  that  the  s-wave  component  of  the  density,  responsible  for 
the  divergent  variance,  gives  no  contribution  to  the  force  and  can  be  therefore  excluded.  The  key 
feature  of  these  estimators  is  their  simplicity;  this  makes  them  straightforward  to  apply  to  many- 
body  systems.  Eor  variational  Monte  Carlo,  this  leads  to  a  direct  way  to  calculate  all  components 
of  ionic  forces  at  once.  In  diffusion  Monte  Carlo  (DMC),  since  the  force  is  calculated  using  the 
electronic  density,  forward  walking  or  reptation  QMC  is  needed  for  unbiased  results.  A  nice 
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feature  of  one  of  our  solutions  is  that  it  retains  the  loeal  eharacter  of  the  original  estimator 
making  its  DMC  implementation  free  of  uncontrolled  approximations.  Our  results  for  various 
molecules  (CH4,  NH3,  H2O,  Li2  H2)  are  more  accurate  than  any  other  method  and  are  the  first 
QMC  results  for  polyatomic  molecules.  We  have  also  worked  on  integrating  a  DFT  code  with 
our  QMC  code  in  order  to  make  the  search  for  stable  structures  as  automatic  as  possible.  Nuclear 
configurations  are  moved  according  to  the  stochastic  gradient  approximation.  The  results  of  this 
work  were  presented  at  meetings  and  a  paper  has  recently  been  submitted. 

Applications 

An  important  part  of  our  work  is  application  to  real  problems  of  current  interest.  The 
first  step  in  this  work  is  to  carry  out  DFT  calculations;  which  can  often  provide  accurate 
predictions  and  can  be  used  to  generating  trial  functions  (and  trial  forces)  in  QMC  calculations. 
During  this  period  (in  work  partially  supported  by  another  grant)  we  have  brought  to  fruition  our 
work  on  nitrogen  at  high  pressure  and  temperature,  as  well  as  predictions  of  new  structures  of 
nitrogen  that  are  metastable  at  low  pressures  and  temperatures.  Nitrogen  is  a  component  of 
many  high-energy  density  materials,  and  elemental  nitrogen  itself  may  be  meta-stable  in  high- 
energy  phases.  Nitrogen  at  high  pressure  has  presented  a  major  challenge  for  many  years,  with 
large  disagreement  of  theory  and  experiment  [5].  However,  recent  experiments[6]  have  observed 
transitions  and  suggest  the  exciting  possibility  that  it  is  sufficiently  meta-stable  that  it  could  be  a 
candidate  energetic  material  with  very  high  energy  per  unit  mass.  Our  simulations  find 
competing  meta-stable  structures  separated  by  large  energy  barriers,  consistent  with  experiment. 
The  most  exciting  result  is  described  in  a  paper  that  has  been  submitted  [7]  -  a  hexagonal 
metallic  chain  structure  that  is  readily  derived  from  the  known  high-pressure  molecular  s  phase 
and  energetically  equivalent  to  the  lowest  energy  high-pressure  structure  found  before. 
Comparison  of  the  energy  with  other  structures  is  shown  in  Fig.  1. 


Figure  1.  Energy  vs.  volume  for  the  known  molecular  phases  of  nitrogen,  the 
previously  predicted  cubic  gauche  phase,  and  the  new  phases  found  in  our  work. 
The  hexagonal  chain  is  metallic  and  is  essentially  the  same  or  more  stable  than 
previously  studied  phases. 
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Our  work  resolves  a  major  puzzle  and  shows  the  importance  of  performing  signification 
ab  initio  simulations  that  revealed  the  new  structure.  These  were  not  observed  in  work  that  was 
limited  to  structures  derived  from  a  simple  cubic  form.  This  work  was  presented  at  APS  March 
(2004)  [8].  These  simulations  were  done  with  the  SIESTA  code[9]  which  is  an  efficient  local 
orbital  code  capable  of  accurate  density  functional  simulations.  This  is  important  because  this  is 
perhaps  the  most  appropriate  code  available  at  present  to  treat  the  condensed  forms  of  complex 
energetic  materials. 

We  have  also  carried  out  studies  of  oxygen  in  molecular  O2  and  ozone  O3  forms,  and 
various  non-molecular  crystalline  forms.  The  former  calculations  are  meant  to  better  understand 
the  structure  of  oxygen  and  ozone  in  molecular  solid  forms,  and  the  test  the  methods  using  the 
SIESTA  code  applied  to  molecular  crystals.  The  simulations  of  high-pressure  non-molecular 
phases  such  as  bee,  fee,  and  sc  have  been  carried  out  to  provide  input  that  can  be  used  to 
determine  universal  force  field  models.  During  the  summer  of  2004  we  will  results  [10]  to  Prof. 
Brenner  and  we  anticipate  continued  interactions  with  MURI  participants  Ammon,  Brenner,  and 
Thompson.  One  interesting  result  is  the  large  effects  of  magnetic  transitions  that  lead  to  large 
changes  in  energy  that  modify  the  order  of  structures  at  intermediate  pressure  ranges.  This 
presents  difficulties  for  the  empirical  models  that  are  yet  to  be  resolved.  They  calculations  are  in 
fact  relevant  to  oxygen  at  very  high  pressure,  and  they  show  that  simple  cubic  (sc)  is  the  favored 
structure  for  a  large  range  of  pressures. 

We  have  started  work  on  molecular  crystals  of  high  energy  density  materials,  with  the 
goal  to  correlate  the  electronic  properties  of  the  molecular  crystal  -  in  particular,  the  band  gaps  - 
with  to  impact  sensitivity  of  nitro  aromatic  explosive  materials. 

Future  Work 

We  will  complete  our  on-going  investigation  of  nitrogen  at  high  pressure  and 
temperature.  Nitrogen  is  a  component  of  many  high-energy  density  materials,  and  elemental 
nitrogen  itself  may  be  meta-stable  in  high-energy  phases.  Our  simulations  find  competing  meta¬ 
stable  structures  separated  by  large  energy  barriers.  The  most  exciting  result  is  a  new  hexagonal 
metallic  chain  structure  that  is  readily  derived  from  the  known  high-pressure  molecular  8  phase 
and  energetically  equivalent  to  the  lowest  energy  high-pressure  structure  found  before.  This 
explains  many  aspects  of  recent  experiments,  which  indicate  metastable  structures,  and  opens  the 
possibility  of  metallic  nitrogen.  A  paper  has  been  submitted  and  we  will  finalize  the  publication. 

We  are  now  carrying  out  full  QMC  calculations  of  predicted  structures  of  nitrogen  to 
verify  if  they  are  indeed  stable.  This  work  is  in  progress  and  detailed  comparisons  with  the  DET 
predictions  will  establish  the  most  accurate  results  possible  by  any  theoretical  method  available 
today.  Preliminary  work  reported  at  the  MURI  meeting  in  October,  2003,  shows  that  the 
structures  are  in  general  agreement  with  the  DET  predictions;  however,  we  have  not  yet  reached 
the  final  conclusions  due  to  difficulties  in  converging  the  energies  with  the  extreme  precision 
needed  to  quantify  small  energy  differences  between  very  forms  of  nitrogen  with  completely 
different  bonding  and  electronic  band  structure. 

We  will  continue  work  on  oxygen  in  molecular  O2  and  ozone  O3  forms,  and  on 
crystalline  forms  to  provide  further  input  for  universal  force  field  models,  in  conjunction  with  the 
work  of  Professors  Ammon,  Brenner,  and  Thompson.  [10] 

We  will  also  complete  work  started  on  high  energy  density  materials,  to  correlate  both 
structure  and  electronic  properties  of  the  molecular  crystal  to  impact  sensitivity  and  other 
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properties  of  nitro  aromatic  explosive  materials.  This  work  is  in  perogress  and  a  paper  is  being 
written.  [1 1] 

We  will  continue  our  work  to  develop  coupled  electronic-ionic  Monte  Carlo  (CEIMC) 
simulations  [See  D.  Ceperley,  M.  Dewing  and  C.  Pierleoni,  (physics/0207006),  in  proceedings  of 
SIMU  conference,  2001,  Springer- Verlag,  eds.  Nielaba  et  al.].  In  this  method,  one  moves  the 
ions  classically  while  calculating  the  electronic  potential  (i.e.  Born-Oppenheimer)  with  QMC. 
The  outstanding  problem  in  this  and  other  Monte  Carlo  methods  is  the  calculation  of  forces.  We 
are  developing  a  new  approach  for  force  estimators  based  on  the  observation  that  the  s-wave 
component  of  the  density,  responsible  for  the  divergent  variance,  gives  no  contribution  to  the 
force  and  can  be  therefore  excluded.  The  key  feature  of  these  estimators  is  their  simplicity;  this 
makes  them  straightforward  to  apply  to  many-body  systems.  We  have  also  worked  on  integrating 
a  DFT  code  with  our  QMC  code  in  order  to  make  the  search  for  stable  structures  as  automatic  as 
possible.  Currently  we  are  writing  up  results  on  nitromethane  (CH3ONO)  which  will  continue  in 
this  period.[12]  We  will  work  toward  combined  DFT/QMC  simulations  of  such  materials. [13] 
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Ab-initio  Potential  Energy  Surfaces  and  Gradients  for  Chemical  Reactions 

(Rodney  J.  Bartlett,  Tom  Hughes,  Euis  Galiano,  and  Joshua  McClellen 
University  of  Elorida,  Quantum  Theory  Project,  Gainesville,  Elorida) 

Goals: 

To  develop  and  apply  quantum  chemical  methods  to  provide  reliable  information  about 
potential  energy  surfaces  for  energetic  molecules,  undergoing  decomposition  and  reaction.  This 
includes  detailed  studies  of  unimolecular  decomposition  and  the  alternative  paths  that  will  likely 
be  encountered  with  additional  interactions  in  the  condensed  phase.  Systems  of  interest  include 
nitramine,  dimethylnitramine,  nitromethane,  methylnitrite,  RDX,  TNAZ,  and  Pox  7.  The  hope  is 
to  provide  accurate  results  for  potential  energy  surfaces  that  can  be  used  to  define  classical 
potentials,  primarily  by  Don  Brenner,  and  be  used  in  simulations  by  Don  Thompson  and  Betsy 
Rice. 

To  accomplish  the  above,  we  use  coupled-cluster  methods  augmented  by  new 
developments  we  have  made  for  larger  molecules,  like  our  new  singular  value  decomposition 
(SVD)  approach  that  leads  to  ‘compressed’  coupled-cluster  equations;  and  the  use  of  frozen 
natural  orbitals  (NO)  to  reduce  the  effective  basis  set  for  large  scale  CC  calculations.  The  NO- 
CC  approach  requires  the  development  of  analytically  computed  forces  (gradients)  to  enable  the 
effective  search  of  potential  energy  surfaces  for  structure  and  transition  states.  We  are  also 
exploring  new  density  functional  methods  to  assess  their  ability  to  describe  energetic  systems,  at 
a  lower  cost. 

We  are  also  developing  a  ‘transfer  Hamiltonian’  that  is  meant  to  be  a  low-rank  operator 
that  enables  the  accuracy  of  CC  theory  to  be  maintained,  but  for  much  more  complicated 
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systems.  In  this  way,  we  hope  to  be  able  to  do  simulations  with  quantum  mechanical  forces  that 
can  be  generated  ‘on-the-fly.’ 

Progress 

We  have  done  a  number  of  calculations  on  nitramine  and  dimethylnitramine  and  related 
systems  to  assess  the  errors  in  the  calculations  and  explore  basis  set  extrapolation  procedures, 
particularly  for  forces.  Many  energy  extrapolation  methods  have  been  used,  but  for  entirely 
satisfactory  results  it  is  equally  important  to  know  the  forces  at  the  basis  set  limit  to  converge 
transition  state  structures  and  associated  activation  barriers. 

We  introduced  our  singular  value  decomposition  (SVD),  ‘compressed”  coupled-cluster 
method  in  two  papers  (1,2)  where  the  detailed  results  can  be  found.  The  basis  idea  is  that  there  is 
substantial  basis  set  linear  dependency  in  high-level,  correlated  quantum  chemical  calculations. 
SVD  is  a  way  to  remove  much  of  this  dependency  by  identifying  an  alternative  set  of  terms  that 
are  weighted  by  their  numerical  importance.  In  normal  CC  theory,  the  number  of  T2  amplitudes 
IS  0  V  for  i,j,„occupied  and  a,b,...  unoccupied  orbitals.  For  triple  excitations  we  require  0  v  , 
etc.  SVD  leads  to  an  alternative  expansion  in  terms  of  amplitudes  t2(X,Y)  instead  of  t2(ia,  jb), 
where  the  X  index  is  a  weighted  (ai)  and  Y  a  weighted  (bj),  so  that  there  a  many  fewer  indices 
than  in  normal  coupled-cluster  theory.  This  is  particularly  important  as  we  go  to  higher 

o 

amplitudes  like  T3,  where  a  computational  ~n  procedure  can  be  dramatically  reduced. 

Similarly,  the  use  of  frozen  natural  orbitals  can  have  a  similar  effect,  but  to  use  them 
effectively,  requires  the  development  of  new  analytical  gradient  methods  for  such  calculations  to 
permit  structures  and  transition  states  to  be  obtained.  We  have  recently  worked  on  this 
problem(3). 

The  concept  of  a  ‘transfer  Hamiltonian’  has  been  developed,  formally  from  wavefunction 
theory,  by  showing  that  an  exact  one -particle  operator  exists  that  can  provide  the  exact  ionization 
potentials,  exact  energy,  and  exact  density  matrix  for  a  molecule  (4).  In  the  case  of  exact  dft,  by 
using  the  density  as  the  critical  quantity,  we  know  that  the  one-particle  operator,  the  Sham 
Hamiltonian,  will  provide  the  molecule’s  exact  density,  its  ionization  potential  as  the  negative 
eigenvalue  of  the  highest  occupied  molecular  orbital,  and  energy  from  inserting  the  density  into 
the  energy  functional.  Consequently,  if  we  can  arrange  to  build  such  a  simplified  one -particle, 
typically,  parameterized  Hamiltonian,  by,  eg,  insisting  that  it  reproduce  the  correct  structures, 
densities,  and  bond-breaking  results  in  coupled-cluster  calculations,  and  demonstrate  that  its 
form  saturates  for  small  clusters,  we  then  have  a  tool  that  permits  the  very  rapid  generation  of 
quantum  mechanical  forces  to  drive  molecular  dynamics.  In  particular,  we  are  doing  this  for 
nitromethane,  where  we  perform  CC  calculations  on  the  dissociation  paths,  and  then  go  to  two 
three,...  nitromethanes  until  all  the  parameters  in  the  chosen  form  of  the  Hamiltonian  are  fixed 
regardless  of  the  number  of  nitromethane  molecules;  then  we  should  expect  that  the  Hamiltonan 
so  defined,  would  be  able  to  describe  many  nitromethane  molecules  with  CC  accuracy.  Only 
such  a  tool  can  be  expected  to  provide  reliable,  quantum  mechanical  results  in  a  time-frame 
suited  to  MD.  The  seamless  transition  between  the  transfer  Hamiltonian  and  the  classical 
potentials  is  particularly  important  for  realistic  sized  condensed  phase  systems. 

Publications: 

•  T.  Kinoshita,  O.  Hino  and  R.J.  Bartlett,  “Singular  value  decomposition 
approach  for  approximate  coupled  cluster  method,”  J.  Chem.  Phys.  119,  7756- 
7762  (2003). 
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•  O.  Hino,  T.  Kinoshita  and  RJ.  Bartlett,  “Singular  value  deeomposition  applied  to  the 
eompression  of  T3  amplitude  for  the  eoupled  eluster  equations,”  J.  Chem.  Phys.  121,  1206-1213 
(2004). 

•  A.  Taube,  A.  Perera  and  R.  J.  Bartlett,  to  be  published. 

•  A.  Beste  and  R.  J.  Bartlett,  “Correlateed  Independent  Partiele  Method;  Numerieal  Results,”  J. 
Chem.  Phys.,  submitted 

Presentations: 

•  Sept  2004  -  Eleetronie  Strueture:  Prineiples  and  Applieations  (ESPA-2004),  Valladolid,  Spain 

•  May  2004  -  Response  Theory  and  Moleeular  Properties,  Sandbjerg  Manor,  Sonderborg, 
Denmark 

•  Eeb  2004  -  “Correlated  Orbital  Theories:  On  wavefunetion  theory  and  density  funetion  theory, 
the  development  of  an  exaet  one-partiele  theory  for  moleoules,”,Theory  and  Applieations  of 
Computational  Chemistry,  Gyeongji,  Korea 

Deeember  2003  -  “Atomie  Seale  Materials  Design:  Modeling  &  Simulation”  Materials  Researeh 
soeiety,  Boston,  MA 

•  Oetober  2003  -  ‘Mb  initio  Predietions  of  PES  for  Chemieal  Reaetions,”  Review  of  Energetie 
Materials  DURINT  and  MURI  Programs,  Aberdeen,  MD 

•  September  2003  -  “Metastable  Moleeules  in  Ground  and  Exeited  States,”  2"*^  Advaneed 
Energeties  Teehnieal  Exehange,  Aberdeen,  MD 

•  September  2003  -  “Coupled-eluster  Methods  and  Their  Applieations  to  Energetie  Moleeules,” 
226*  ACS  National  Meeting,  New  York,  NY 

•  August  2003  -  “High  level  Couple  Cluster  Theory:  What  Did  We  Learn?”  8*  European 
Conferenee  on  Quantum  Systems  in  Chemistry  and  Physies,  Spetses,  Greeee 

•  July  2003  -  “Hh  Initio  Density  Eunetional  Theory,”  Eleetron  Correlation:  Ab  initio  Methods 
and  Density  Eunetional  Theory,  Satellite  Meeting  of  the  Xlth  International  Congress  of  Quantum 
Chemistry,  Bad  Herrenalb,  Germany 

•  April  2003  -  “Erontiers  in  Theoretieal  Chemistry,”  a  Symposium  in  Honor  of  Prof.  Rudolph 
A.  Mareus,  Los  Angeles,  CA 

Mareh  2003  -  “Erom  Wave  Eunetion  Theory  to  Density-Eunetional  Theory  and  Baek,”  225* 
ACS  National  Meeting,  New  Orleans,  LA 

•  Oetober  2002  -  “Sealable  Software  for  Computational  Chemistry,”  Grid  Computer 
Conferenee,  University  of  Kentueky,  Lexington,  KY 

•  Oetober  2002  -  “Hh  initio  Predietions  of  PES  for  Chemieal  Reaetions,”  MURI  Kiek-off 
Meeting,  Aberdeen,  MD 

•  September  2002  -  “Large  Seale  Dynamies  with  Quantum  Meehanieal  Eorees,”  Symposium 
and  Summer  Sehool  on  Nano  and  Giga  Challenges  in  Mieroeleetronies  Researeh  and 
Opportunities,  Moseow,  Russia 

•  September  2002  -  “Predietive  Theory  from  Moleeules  to  Materials,”  Seienee  at  the  Edge, 
Miehigan  State  University 

•  July  2002  -  “Advanees  in  Eleetronie  Strueture  Theory:  Current  Trends  and  Euture  Prospeets,” 
International  Conferenee  of  Theoretieal  Chemieal  Physies,  IV,  Marly-le-Roi,  Eranee 

•  June  2002  -  “Predietive  Theory  from  Moleeules  to  Materials,”  Symposium  to  Initiate  Joint 
Ph.D.  Program  between  Eranee  and  the  United  States,  Strasbourg,  Eranee 

•  June  2002  -  “Large  Seale  Simulations  with  Quantum-Meehanieal  Eorees,”  European  Materials 
Researeh  Soeiety  Spring  Meeting,  Strasbourg,  Eranee 


14 


Students  Supported  by  MURI: 

•  My  most  recent  Ph.D.  DeCarlos  Taylor,  is  now  on  an  NRC  postdoctoral  position  at  the  Army 
Research  Laboratory. 

•  Several  other  students  have  been  partly  supported  by  this  MURI.  Dr.  Ariana  Beste,  Ph.D., 
2004,  now  at  Oak  Ridge  national  Laboratory;  Dr.  Tom  Henderson,  Ph.D.,  2004,  now  at  the 
NMRC,  Cork,  Ireland;  Dr.  Anthony,  Yau,  Ph.D.,  2004,  now  at  ACES  QC,  Gainesville,  FL,  and 
current  students:  Josh  McClellan,  Tom  Hughes,  and  Luis  Galiano. 

Interacgtions  with  Army  Scientists: 

Dr.  DeCarlos  Taylor 
Dr.  Steve  Bunte 
Dr.  Betsy  Rice 
Dr.  Shashi  Kama 
Dr.  Cary  Chabalowski 


Quantum-Based,  Reactive  Potentials  for  Simulating  Shock  Dynamics  of  Condensed-Phase 
Energetic  Materials 

(Donald  W.  Brenner,  Materials  Science  and  Engineering  Dept.,  North  Carolina  State  University, 
Raleigh,  N.C.) 

We  are  currently  developing  a  transferable  and  robust  analytic  reactive  potential  for  C,  H, 
O  and  N  containing  species  that  is  valid  in  both  the  molecular  and  condensed  phase  regimes  and 
that  allows  chemical  reactivity  through  all  changes  of  detonation  and  shock  propagation,  as  well 
as  slower  deflagration  processes.  The  analytic  functional  form  is  based  on  a  chemically  sound 
bond-order  formalism  that  has  proven  extremely  powerful  for  describing  reactivity  in 
hydrocarbon  systems.  Accurate  first  principles  quantum-mechanical  calculations  are  being  used 
to  both  determine  appropriate  functional  forms  and  parameters  entering  this  formalism,  and  to 
validate  specific  chemical  reaction  paths  and  rates  for  unimolecular  dissociation  and 
recombination  produced  by  the  analytic  potential.  These  computationally  efficient  potentials  are 
enabling  large-scale,  three-dimensional  molecular  dynamics  simulations  to  be  carried  out  at  the 
Army  Research  Eaboratory  that  will  predict  system  properties  related  to  shock  initiation  and 
detonation  of  a  wide  range  of  both  existing  and  potentially  safer  and  more  powerful  high 
explosives  (i.e.  molecular  design). 

Progress: 

A  reactive  empirical  bond  order  potential  for  NxHy  species  has  been  completed,  and  the 
functions  and  parameters  sent  to  Dr.  Betsy  Rice  at  the  ARE  for  incorporation  into  the  current 
DoD  molecular  shock  modeling  codes.  Dr.  Rice  has  completed  the  initial  coding  needed  to 
derive  forces  from  the  potential  energy  expression,  and  will  initiate  molecular  modeling 
simulations  with  this  code  and  forces.  This  potential  will  initially  be  used  in  conjunction  with  Dr. 
Rice  and  Prof.  Thompson  to  model  shock  chemistry  of  liquid  hydrazine. 

This  version  of  the  NxHy  potential  is  consistent  in  form  with  our  CxHy  potential 
previously  developed,  and  consistent  with  a  new  potential  for  CHO  species  developed  by  Prof 
Sinnott  at  the  University  of  Florida.  We  have  begun  building  on  this  common  functional  form  by 
combining  it  with  new  bonding  and  charge  screening  screening  approaches  developed  as  part  of 
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our  MURI  effort  (see  prior  report).  Currently,  effeetive  pair  terms  in  this  many-body  potential 
have  been  developed  for  all  pairs  of  atom  types  exeept  for  N-0. 

The  fitting  data  base  for  our  NxHy  potential  is  extensive,  and  ineludes  bond  energies  and 
bond  lengths  for  the  moleeular  speeies  H2,  N2,  NH,  NH2,  NH3,  HNNH,  H2NNH2,  H3NNH3, 
HNNH3,  the  energy  barrier  for  the  inversion  of  NH3,  and  eondensed  phase  properties  of  liquid 
ammonia.  This  is  the  first  potential  energy  expression  to  address  the  properties  of  eaeh  of  the 
speeies  in  a  single  funetional  form,  a  feature  that  will  be  key  in  utilizing  the  potential  to 
eharaeterize  the  properties  of  hydrazine.  Beeause  many  of  the  properties  of  hydrazine  are 
reasonably  well  known  from  experiment,  the  results  of  simulations  using  our  potential  will  both 
help  validate  the  potential,  as  well  as  lend  important  new  insights  into  the  reaetive  dynamies  of 
these  systems. 

The  Sinnott  oxygen  potential  was  not  fit  to  nor  validated  against  high-pressure  phases  of 
oxygen  that  may  play  a  role  as  transient  intermediates  in  detonating  solids.  To  validate  the 
potential  for  this  regime.  Prof.  Martin’s  group  has  used  first  prineiples  methods  to  ealeulate  the 
energy  of  oxygen  in  faee  eentered  eubie,  simple  eubie  and  diamond  eubie  lattiees  as  a  funetion  of 
lattiee  eonstants.  Both  spin  restrieted  and  spin  unrestrieted  ealeulations  were  earried  out. 
Comparisons  between  these  ealeulations  and  the  energy-lattiee  eonstant  relations  of  the  Sinott 
potential  are  eurrently  being  quantified. 

An  informal  workshop  was  held  at  North  Carolina  State  University  on  June  29  and  30, 
2004.  Hosted  by  Prof.  Brenner  and  the  NC  State  Department  of  Materials  Seienee  and 
Engineering,  a  purpose  of  this  workshop  was  to  diseuss  eommon  issues  and  progress  related  to 
the  development  of  many-body  interatomie  potential  energy  funetions  by  some  of  the  leading 
university  researehers  in  the  southeast.  The  workshop  was  attended  by  Prof  Sinnott  and  Prof. 
Simon  Philpot  of  the  University  of  Florida,  Prof.  Steve  Stuart  of  Clemson  University,  and  Prof. 
Judith  Harrison  of  the  Naval  Aeademy,  as  well  as  postdoes  and  students  from  these  groups. 
Progress  on  the  hydrazine  potential  being  developed  in  the  MURI  was  diseussed  by  Prof. 
Brenner,  ineluding  the  new  line-of-state  sereening  and  topology  funetions  being  used  to  deseribe 
bonding  range  and  bonding  type  between  pairs  of  atoms.  Prof  Sinnott  diseussed  her  new  CHO 
potential  mentioned  above.  Prof.  Stuart  deseribed  his  reeent  work  on  variable  eharge  methods  in 
whieh  an  effeetive  eharge  is  assigned  to  an  atom  aeeording  to  its  environment  using  constraints 
related  to  the  relative  electronegativity  and  hardness  of  an  atom,  and  how  these  methods  can  be 
incorporated  into  bond  order  potentials.  Prof.  Simon  Phillpot  discussed  his  numerical  method  for 
calculating  long-range  Coulomb  interactions  in  which  a  convenient  cut-off  radius  is  chosen  and 
any  residual  charge  is  placed  at  a  sphere  with  the  cut-off  radius.  The  method  is  accurate  and 
convenient  to  code  as  an  independent  subroutine.  We  are  currently  exploring  the  utility  methods, 
including  their  incorporation  into  the  ARL  modeling  codes. 

Supported  Personnel: 

Ben  Gilbert,  graduate  student 

Kai  Wang,  research  assistant  professor 

Beyond  Steady-State  Deflagration  Models:  Molecular  Dynamics  at  a  Gas-Liquid  Interface 
(John  E,  Adams,  T.  Szabo,  and  A,  Siavosh-Haghighi;  University  of  Missouri,  Columbia, 

MO  65211) 

State-of-the-art  models  of  the  steady-state  deflagration  of  an  energetic  material'  have 
proved  reasonably  successful  in  reproducing  experimental  combustion  rates  in  a  number  of 
systems.  Inasmuch  as  these  models  are  based  on  continuum  descriptions  of  the  condensed 
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phases,  they  are  likely  to  have  reached  the  limit  to  which  they  can  be  extended  and  enhanced, 
however,  in  the  absence  of  detailed  molecular  information.  Furthermore,  they  necessarily  are 
applicable  only  to  systems  in  which  a  steady  state  has  been  established  and  thus  are  of  little 
direct  help  in  modeling  the  initiation  of  combustion.  Understanding  how  a  combusting  system 
reaches  steady  state  will  require  additional  insight  into  the  mass  and  energy  balance  during  the 
early  stages  of  combustion,  with  part  of  that  picture  being  the  transfer  of  energy  from  nascent 
gaseous  combustion  products  to  the  thin  liquid  film  that  forms  at  the  surface  of  a  burning  solid. 
Even  within  the  context  of  the  extant  steady-state  models,  though,  there  is  a  need  for  molecular- 
level  data  for  the  implementation  and  detailed  validation  of  these  models,  data  such  as 
temperature-dependent  vapor  pressures  and  a  characterization  of  the  gaseous  product  dynamics 
in  the  presence  of  temperature  gradients. 

Our  initial  focus  has  been  on  obtaining  a  basic  understanding  of  the  energy  transfer  from 
translationally  hot  gas  atoms  to  a  prototypical  liquid  surface.  The  simplest  non-trivial  system 
that  nonetheless  has  the  potential  of  capturing  the  gross  dynamics  of  a  more  complex  system 
consists  of  atoms  interacting  through  Lennard- Jones  forces,  a  system  such  as  the  one  investigated 
previously  by  Tribe  et  al.  Simulations  of  this  system  have  been  performed  using  the  modular 
DL_POLY_2  code^,  with  additional  modules  being  developed  as  necessary  to  handle  system 
initialization  and  the  final  energy  analysis. 

We  begin  by  equilibrating  a  liquid  sample  at  a  selected  temperature  {NVT  ensemble, 
periodic  boundary  conditions,  cubic  simulation  cell)  then  create  an  interface  (actually,  two 
interfaces)  by  expanding  the  simulation  cell  length  in  the  z  dimension.  (Periodic  boundary 
conditions  in  three  dimensions  are  maintained  following  the  expansion.)  At  this  point,  we  can 
either  (1)  proceed  at  constant  energy — the  system  then  cools  as  kinetic  energy  is  converted  to 
potential  energy — ^until  a  final  average  temperature  is  established  or  (2)  continue  the  simulation 
at  the  desired  final  temperature.  (Prior  to  this  final  equilibration,  the  net  center-of-mass  velocity 
of  the  liquid  sample,  which  is  an  artifact  of  the  way  in  which  the  interfaces  are  created,  must  be 
substracted.)  These  procedures  lead  to  a  density  contour  in  the  system  such  as  that  shown  here. 


Note  that  as  the  temperature  of  the  sample  increases,  the  width  of  the  interfaces  also  increases. 
Furthermore,  at  the  highest  temperature  shown,  one  easily  sees  a  finite,  uniform  vapor  density  in 
the  volume  of  the  simulation  cell  that  initially  was  empty.  The  initial  configurations  for  the 
energy  transfer  studies  are  then  generated  by  recording  snapshots  of  this  system  periodically  as 
the  NVT  simulation  proceeds. 
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For  the  actual  studies  of  energy  transfer,  we  direct  atoms  with  a  well-defined  kinetic 
energy  at  the  surface  configurations  generated  as  described  above,  the  initial  positions  of  these 
atoms  being  chosen  at  random  within  a  plane  parallel  to  the  liquid  surface  and  at  a  distance 
sufficiently  far  from  it  that  there  is  no  initial  interaction.  The  energies  (and  scattering  angles)  of 
these  atoms  are  then  analyzed  at  the  point  that  (1)  they  rebound  from  the  surface  and  reach  the 
plane  from  which  they  initially  were  injected  into  the  system  (“scattered  atoms”)  or  (2)  the 
trajectory  times  exceed  a  maximum  value,  beyond  which  the  atoms  are  assumed  to  be  trapped 
and  equilibrated  with  the  surface  (“trapped  atoms”).  However,  if  at  the  maximum  trajectory  time 
an  atom  is  moving  away  from  the  surface  and  experiencing  no  interactions  with  other  species  in 
the  system,  then  the  atom  is  also  counted  being  “scattered”. 

The  final  energy  distribution  of  initially  hot  atoms  directed  towards  a  liquid  surface  (at  a 
Lennard-Jones  reduced  temperature  of  0.719;  the  incident  angle  is  55°  with  respect  to  the  surface 
normal)  is  given  in  the  figure  below.  Note  that  the  full  distribution  (labeled  “All”)  is  strongly 
peaked  at  low  energies,  a  result  that  reflects  the  thermal  equilibration  of  trapped  species  with  the 
surface.  (The  average  kinetic  energy  of  atoms  equilibrated  with  the  surface  for  this  system 
would  be  5.44  kJ/mol.) 


■ 


Final  Energy  (kJ/mol) 


It  is  interesting  to  compare  these  results  with  those  reported  by  Nathanson  and  co- 
workers  ,  who  investigated  atom  scattering  from  a  liquid  metal  surface.  (The  parameters  in  the 
present  study  have  been  chosen  to  mimic  those  experiments.)  In  practice,  of  course,  energy 
analysis  in  the  experiments  necessitates  the  use  of  a  detector  having  a  finite  entrance  aperture. 
Thus,  if  one  models  the  experiment  by  analyzing  the  final  energies  of  atoms  scattered  at  the 
specular  angle  (and  using  a  relatively  large  entrance  aperture  to  improve  the  counting  statistics  ), 
one  finds  a  very  different  energy  distribution,  one  that  is  peaked  at  about  50%  energy  transfer, 
with  an  average  energy  transfer  that  is  somewhat  less  than  that. 
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The  key  to  understanding  the  relationship  between  collisional  energy  transfer  and  the 
establishment  of  a  steady-state  temperature  in  a  deflagrating  energetic  material  is  the  dependence 
of  the  energy  transfer  on  the  surface  temperature.  The  figure  below  shows  a  comparison  of  the 
scattered  energy  distributions  calculated  at  two  different  temperatures  (reduced  temperatures  of 
0.719  and  0.923).  The  distributions  in  fact  are  found  to  be  quite  similar,  but  it  is  important  to 


note  that  these  distributions  have  been  scaled  in  each  case  by  the  number  of  scattered  atoms.  The 
significant  difference  between  the  two  sets  of  results  is  that  at  the  lower  temperature  80%  of  the 
atoms  are  scattered  rather  than  trapped,  while  at  the  higher  temperature  this  percentage  drops  to 
64%.  This  increased  incidence  of  trapping  (and  thus  thermal  equilibration)  as  the  surface 
temperature  increases  can  be  correlated  with  greater  surface  roughness — the  broader  surface 
contour  at  the  higher  temperature  also  testifies  to  the  temperature-dependence  of  the  surface 
topology.  (Surface  roughening  was  also  posited  by  King  et  ah,  to  explain  the  experimental 
scattering  of  argon  from  perfiuorinated  polyethers."^)  A  complete  analysis  of  the  surface 
temperature  dependence  of  energy  transfer  in  our  prototype  system  is  in  its  final  stages,  with  a 
manuscript  describing  the  work  being  prepared  for  submission  by  the  end  of  the  year. 

We  also  are  pursuing  an  investigation  of  a  practical  energetic  material,  nitromethane,  for 
which  potential  parameters  and  a  characterization  of  the  melting  behavior  are  available  in  the 
work  of  Agrawal,  Rice,  and  Thompson.^  These  calculations  are  significantly  more  demanding 
than  are  those  relevant  to  the  Lennard-Jones  system  in  that  we  include  all  internal  degrees  of 
freedom  of  the  nitromethane  molecules.  (In  an  extension  of  the  present  work,  we  will  freeze  the 
internal  coordinates  of  the  molecules  and  assess  the  extent  to  which  the  collisional  energy 
transfer  is  altered  when  molecular  vibrations  no  longer  provide  an  accessible  energy  sink.)  The 
surface  contour  in  this  system  is  given  in  the  figure  below. 
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While  this  work  is  still  ongoing,  an  example  of  the  system 
configuration  is  shown  here.  The  important  point  to  recognize  in  examining 
this  configuration  is  that  the  surface  layer  is  by  no  means  simple,  that  it  is 
rough  on  a  microscopic  scale.  Accordingly,  an  atom  colliding  with  this 
surface  is  expected  to  experience  a  high  probability  of  trapping.  (In  an 
example  trajectory  computed  for  this  system,  64%  of  the  incident  energy  was 
transferred  as  a  result  of  the  collision.)  We  expect  to  have  a  manuscript 
describing  our  initial  work  on  this  system,  work  that  will  include  a 
determination  of  the  equilibrium  vapor  pressure  of  nitromethane  above  the 
liquid  surface,  ready  for  submission  by  year  end. 

Finally,  we  have  begun  looking  at  ozone  as  a  system  that  will  be 
useful  in  connecting  the  simulations  with  extant  continuum  models  of 
deflagration.  One  would  like  to  use  a  potential  model  in  this  system  that 
correctly  accounts  for  all  species  present  in  the  system  during  burning. 
Work  leading  to  a  potential  model  with  the  desired  properties  is  presently 
under  development  by  the  Brenner  group.  We  expect  to  have  access  to  these 
results  shortly. 
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Development  of  Methods  for  Predicting  Solvation  and  Separation  for  Energetic 
Materials  in  Supercritical  Fluids 

(Christopher  J.  Cramer,  Jason  Thompson,  Casey  P.  Kelly,  Benjamin  J.  Lyneh  and  Donald  G. 
Truhlar;  Department  of  Chemistry  and  supercomputing  Institute,  University  of  Minnesota 
Minneapolis,  Minnesota) 

The  current  goal  of  the  MURI  work  here  at  the  University  of  Minnesota  is  to  aehieve  a 
better  understanding  of  the  solubility  and  other  properties  of  substances  in  supercritical  fluids 
and  to  employ  that  understanding  in  the  development  of  supereritical  fluid  teehnologies  for 
recycling  and  reclamation  of  energetic  materials.  In  particular,  we  plan  to  model  solubility  in 
supercritical  fluids  because  this  property  has  important  applications  to  high-energy  materials. 
We  are  developing  the  tools  and  software  for  accurate  predictive  methods  and  models  for 
solvation  of  energetic  materials  in  supereritical  carbon  dioxide;  we  are  implementing  them  into 
easy-to-use  eodes  that  are  freely  available  via  the  World  Wide  Web.  Thus  Army  and  MURI 
researehers  will  be  able  to  use  these  models  and  methods  to  facilitate  the  design  of  practical 
procedures  for  extraction,  recycling,  and  reusing  materials. 

Sinee  the  MURI  projeet  began  in  2002,  we  have  focused  on  extending  existing 
eontinuum  solvation  models  to  supercritical  solvents.  Our  existing  eontinuum  solvation  models 
(i.e.,  solvation  models  for  ordinary  liquid  solvents)  utilize  the  generalized  Born  (GB)  model  in  a 
self-consistent  reaction  field  step  to  aecount  for  long-range  electrostatic  effects.  Note  that  the 
GB  model  represents  the  solute  as  a  collection  of  atom-centered  spheres  and  atom-centered  point 
charges.  Our  model  also  aceounts  for  solvation  phenomena  beyond  eleetrostatics,  including 
cavitation,  dispersion,  hydrogen  bonding,  and  struetural  changes  taking  plaee  in  the  first 
solvation  shell(s),  as  well  as  loeal  ehanges  in  the  permittivity  of  the  solvent,  by  modeling  these 
effects  as  being  proportional  to  the  solvent-aecessible  surface  areas  (SASAs)  of  the  atoms  in  the 
solute.  The  eonstants  of  proportionality  are  called  atomic  surface  tension  coefficients.  They 
contain  parameters  (ealled  atomic  surface  tension  parameters)  that  are  optimized  against  a 
training  set  of  experimentally  known  free  energies  of  solvation  in  water  and  organic  solvents  and 
that  depend  on  a  set  of  experimental  solvent  deseriptors  of  the  solvent.  Our  initial  approaeh  for 
developing  a  continuum  solvation  model  for  supercritical  CO2  is  to  treat  the  critieal-fiuid  solvent 
as  a  continuum  liquid,  i.e.,  a  continuous  and  homogenous  medium  eharacterized  by  the  bulk 
dielectric  constant,  augmented  with  surface  tension  coeffieients  optimized  for  liquid  solvents, 
and  then  optimize  a  set  of  solvent  deseriptors  for  supereritical  CO2  that  are  functions  of 
temperature  and  pressure. 

In  the  rest  of  this  doeument,  a  summary  of  the  progress  made  for  this  MURI  project  and  a 
list  of  publications  and  manuscripts  in  preparation  of  MURl-supported  work  are  given. 

Progress: 

Improved  Charge  Models 

During  the  first  year  and  a  half  of  the  MURI  project,  we  foeused  on  developing  new 
eharge  models  that  provide  aecurate  and  reliable  charges  to  use  in  the  GB  model.  Beeause 
diffuse  functions  may  be  important  for  an  even-  handed  description  of  conformational  energies 
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and  of  negative  functional  groups  found  in  many  of  the  compounds  of  interest,  we  have 
developed  a  new  population  analysis  method,  called  redistributed  Lowdin  population  analysis 
(RLPA).  This  method  provides  partial  atomic  charges  that  are  less  sensitive  to  the  inclusion  of 
diffuse  functions  to  the  basis  set  than  Lowdin  population  analysis  partial  atomic  charges.  We 
have  also  created  a  new  class  IV  charge  model,  called  charge  model  3  (CMS).  The  CMS  charge 
model  uses  a  semiempirical  charge-mapping  scheme  that  is  a  function  of  the  Mayer  bond  orders 
in  the  molecule.  This  scheme  systematically  corrects  errors  in  bond  dipoles  calculated  from 
low-level  charges,  such  as  Lowdin  population  analysis  charges.  The  CMS  charge  model  has 
several  improvements  over  previous  class  IV  charge  models  developed  in  our  research  group.  In 
particular,  the  CMS  parameters  were  optimized  against  a  larger  and  more  diverse  training  set 
than  previous  charge  models.  Second,  when  diffuse  basis  functions  are  used,  CMS  maps  RLPA 
charges  rather  than  Lowdin  population  analysis  charges.  In  addition,  the  CMS  parameters  were 
chosen  so  that  the  resulting  charges  are  not  as  sensitive  to  variations  observed  in  the  Mayer  bond 
order  when  diffuse  basis  functions  are  included  in  the  basis  set.  Finally,  for  some  wave 
functions,  CMS  uses  a  new  charge-mapping  scheme  for  solutes  containing  both  N  and  O. 
Validation  of  Charge  Model  3 

Our  progress  in  2003  and  at  the  beginning  of  2004  included  the  validation  of  the  CMS 
charge  model  for  predicting  accurate  charge  distributions  of  high-energy  materials  (HEMs)  and 
compounds  analogous  to  them.  Charge  model  3  is  a  parameterized  model,  where  the  parameters 
have  been  optimized  against  a  training  set  of  dipole  moment  data.  The  training  set  is  large  (398 
data)  and  diverse,  containing  many  different  types  of  functional  groups  that  one  encounters  in 
organic  chemistry.  However,  the  functional  groups  found  in  HEMs  are  either  under-represented 
with  respect  to  other  functional  groups  present  in  the  training  set  (nitro  compounds)  or  not 
represented  at  all  (nitramines,  for  example)  in  the  CMS  training  set.  To  determine  whether  or  not 
the  CMS  parameters  are  applicable  for  predicting  accurate  charge  distributions  of  HEMs,  we 
assembled  a  test  set  of  compounds,  including  nitramide,  dimethylnitramine  (DMNA), 
1,3,3-trinitroazetidine  (TNAZ),  l,3,5-trinitro-5-triazine  (RDX),  and  hexanitrohexaaza- 
isowurtzitane  (HNIW),  and  compared  the  dipole  moments  for  different  conformations  of  these 
compounds  (14  in  all)  calculated  using  CMS  charges  to  high-level  theoretical  density  dipole 
moments.  Note  that  a  density  dipole  moment  is  calculated  from  the  one-electron  density  as  an 
expectation  value  of  the  dipole  moment  operator,  and  density  dipole  calculations  do  not  provide 
the  partial  atomic  charges  needed  for  condensed-phase  modeling.  These  tests  have  shown  that 
partial  atomic  charges  from  CMS  can  be  used  to  predict  dipole  moments  for  these  types  of 
compounds  with  similar  or  better  accuracy  as  for  compounds  in  the  CMS  training  set. 
Eurthermore,  this  good  agreement  is  obtained  even  when  relatively  inexpensive  wave  functions 
are  used  for  the  solute,  which  is  important  for  larger  solutes  like  HNIW.  Eor  the  above  test  set  of 
nitramines,  we  have  also  shown  that  atomic  charges  calculated  with  CMS  yields  polarization  free 
energies  (computed  from  the  GB  model)  for  the  condensed  phase  that  are  less  sensitive  to  the 
level  of  treatment  of  electron  correlation  than  those  given  by  other  methods.  This  invariance  to 
the  level  of  treatment  of  electron  correlation  demonstrated  by  the  CMS  model  will  be  important 
in  future  work,  which  will  focus  on  the  further  development  and  testing  of  different  theoretical 
methods  for  modeling  the  solid-state  condensation  of  HEMs. 

A  New  Continuum  Solvation  Model  Based  on  Charge  Model  3 

We  have  also  determined  a  new  set  of  atomic  surface  tension  parameters  to  be  used  in 
conjunction  with  GB  model  employing  CMS  charges.  These  new  parameters  will  be  used  in  the 
optimization  of  solvent  descriptors  for  supercritical  CO2.  They  were  optimized  against  2237 
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experimental  free  energies  of  solvation  of  solutes  in  water  and  90  organie  solvents  and  79  free 
energies  of  transfer  between  water  and  12  different  organie  solvents  (the  free  energy  of  transfer 
of  a  solute  defines  the  partition  eoeffieient  of  a  solute  between  two  different  solvents).  The 
resulting  model,  ealled  SM5.43R,  ean  prediet  free  energies  of  solvation  of  solutes  in  water 
solvent  and  in  organie  solvents,  provided  the  solvent  deseriptors  for  the  organie  solvent  of 
interest  are  known.  In  addition,  SM54.3R  predicts  aqueous  free  energies  of  solvation  a  factor  of 
one  (ions)  to  two  (neutrals)  more  accurately  than  the  continuum  solvation  models  available  in  the 
Gaussian  electronic  structure  program  and  free  energies  of  organic  solvation  a  factor  of  six  to 
seven  (!)  more  accurately  than  these  models.  This  represents  a  resounding  success  for  the  new 
model. 

In  2003,  we  initially  optimized  SM5.43R  parameters  for  the  mPWlPW91  (also  denoted 
MPW25)  hybrid  density  functional  method  with  the  6-31G(d)  and  6-31+G(d)  basis  sets.  The 
mPWlPW91  combines  Barone  and  Adamo's  modified  version  of  Perdew  and  Wang's  exchange 
functional,  Perdew  and  Wang's  correlation  functional,  and  a  percentage  X  of  exact  Hartree-Fock 
exchange.  In  2004,  we  carried  out  further  SM5.43R  parameter  optimizations  for  the 
MPWX/MIDI! ,  MPWX/MIDI16D,  and  MPWX/6-3I+G(d,p)  combinations  of  electronic  structure 
method  and  basis  set  with  X  =  0,  25,  42.8,  and  60.6,  and  for  MPWX/6-31G(d)  and 
MPWX/6-31+G(d)  with  X  =  0,  42.8,  and  60.6.  For  each  of  the  five  basis  sets,  we  found  no 
significant  loss  in  the  accuracy  of  the  model  when  parameters  averaged  over  the  four  values  of  X 
are  used  instead  of  the  parameters  optimized  for  a  specific  value  of  X.  This  is  a  significant  result 
because  for  a  given  property  of  a  molecule  or  reaction,  it  may  be  useful  to  optimize  the  value  of 
the  fraction  of  Hartree-Fock  exchange  in  MPWX  With  an  optimized  value  of  X  in  hand  for  a 
particular  problem,  it  is  useful  to  have  a  solvation  model  already  parameterized  for  that  value  of 
X 

Implementations  ofRLPA,  CMS,  and SM5,43R 

We  have  implemented  this  new  model  in  three  quantum  chemistry  computer  programs, 
namely,  gamessplus,  hondoplus,  and  SMXGAUSS,  all  of  which  are  freely  available  on  the 
internet  (see  http://comp.chem.umn.edu/truhlar).  We  have  made  various  optimizations  to 
HONDOPLUS  (and  consequently  to  SMXGAUSS,  which  is  based  on  the  electronic  structure  code 
implemented  in  HONDOPLUS)  to  allow  it  to  run  four  times  faster  for  density-functional  theory 
calculations.  In  addition  to  these  optimizations,  we  have  also  increased  the  portability  of  our 
solvation  codes,  allowing  us  to  make  better  use  of  the  computational  resources  available.  These 
increases  to  speed  and  portability  have  decreased  turnaround  time  to  perform  calculations  for  the 
development  of  new  solvation  models. 

Development  of  a  Training  Set  of  Solvation  Data  of  Solutes  in  Supercritical  CO2 

In  addition  to  the  GB  model,  our  continuum  solvation  model  incorporates  short-range 
solvation  effects  with  semiempirical  atomic  surface  tension  terms.  These  semiempirical  atomic 
surface  tension  terms  contain  parameters  that  are  optimized  against  a  training  set  of  known 
experimental  free  energies  of  solvation.  In  order  to  develop  a  continuum  solvation  model  for 
supercritical  carbon  dioxide,  a  training  set  of  experimentally  determined  free  energies  solvation 
of  various  solutes  in  supercritical  carbon  dioxide  is  required.  Therefore  we  carried  out  a 
thorough  search  of  the  literature  on  relevant  data  for  supercritical  carbon  dioxide.  As  a  result  of 
this  search,  of  our  reading  in  this  area,  and  our  discussions  at  the  meeting  in  Aberdeen  (October 
2002),  we  came  to  the  conclusion  that  our  modeling  input  (for  deriving  new  surface  tension 
parameters)  for  the  supercritical  aspects  of  our  solvation  model  will  have  to  be  based  on 
solubilities  and  vapor  pressures  of  pure  substances  rather  than  free  energies  of  solvation  and 
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partition  coefficients.  In  eontrast,  all  of  our  past  work  was  based  on  starting  with  free  energies  of 
solvation  and  partition  eoeffieients.  So  we  investigated  the  effeet  of  this  on  our  studies  and  our 
modeling  plan. 

First  we  derived  the  equations  relating  solubilities,  vapor  pressures  of  pure  substanees, 
and  free  energies  of  solvation  to  one  another  in  the  ease  of  ideal  solutions  (i.e.,  when  all  aetivity 
eoeffieients  and  fugaeity  eoeffieients  are  unity).  Then  we  ereated  a  test  set  of  85  solutes  (both 
liquids  and  solids,  mainly  eompounds  eomposed  of  H,  C,  N,  and  O,  but  also  a  few  halogenated 
eompounds  to  add  diversity)  for  testing  how  well  the  ideal  solution  equations  hold  for  saturated 
solutions.  This  is  a  very  fundamental  question  in  ehemieal  thermodynamics,  and  this 
investigation  was  essential  to  the  design  of  subsequent  steps  in  our  modeling  effort.  In  addition, 
we  were  unaware  of  any  tests  eomparable  to  our  own  that  have  been  published  in  the  literature. 

For  liquid  solutes,  we  found  that  we  ean  predict  solubilities  by  our  methods  with 
eomparable  aecuraey  to  what  we  have  previously  eome  to  expeet  for  free  energies  of  solvation 
and  partition  eoeffieients.  We  also  found  that  we  ean  use  experimental  vapor  pressures  as  part  of 
our  model  when  available,  but  we  can  equally  well  prediet  the  vapor  pressures  that  are  needed  as 
part  of  the  solubility  ealeulation,  and  it  did  not  degrade  the  aeeuracy  of  our  predieted  solubilities. 
In  addition,  our  predietions  were  almost  as  aeeurate  as  ean  be  done  when  experimental  free 
energies  of  solvation  are  available. 

The  results  for  solid  solutes  were  not  as  aeeurate  as  our  results  for  liquid  solutes,  but  they 
were  still  eneouraging.  Nevertheless,  by  ineluding  both  liquid  and  solid  solutes  in  our  tests,  we 
aehieved  a  better  fundamental  understanding  of  the  issues  involved  in  modeling  solubilities,  and 
we  opened  ourselves  to  the  possibility  of  ineorporating  a  broader  database  for  developing 
parameters,  a  possibility  that  we  will  take  advantage  of  and  that  will  be  a  eritieal  aspeet  of  our 
model  development. 

We  have  started  to  assemble  a  training  set  of  solutes  with  known  experimental 
solubilities  in  supereritieal  CO2  over  a  wide  temperature  and  pressure  range.  In  order  to  ereate  a 
set  of  generally  useful  solvent  deseriptors,  this  training  set  ineludes  both  high-energy  eompounds 
and  eompounds  eontaining  a  variety  of  other  funetional  groups.  In  addition  to  searehing  the 
ehemieal  literature  and  submitting  data  requests  to  the  Army,  we  have  been  in  eontaet  with 
Victor  Stepanov  at  Pieatinny  Arsenal  in  New  Jersey  and  Lev  Krasnoperov  (through  Dr. 
Stepanov)  at  the  New  Jersey  Institute  of  Teehnology  about  measuring  solubilities  of  nitramine 
eompounds  using  in  situ  UV  spectroseopy.  In  addition  we  suecessfully  obtained  institutional 
authorizations  to  obtain  Militarily  Critieal  Teehnieal  Data  (CPIA  M3  and  M4  Manual  Units),  and 
we  purehased  a  safe  to  hold  this. 

Current  Research 

In  order  to  predict  solubilities  of  high-energy  materials  in  supereritieal  CO2,  we  will  need 
to  be  able  to  predict  the  pure-solute  vapor  pressure  of  these  eompounds.  The  required  solvent 
deseriptors  for  these  types  of  eompounds  are  not  experimentally  known,  and  so  we  need  a 
method  to  estimate  them.  Platts  and  eoworkers  have  developed  a  model  that  ean  be  used  to 
estimate  two  (Abraham's  aeidity  and  basieity  parameter)  of  the  four  solvent  deseriptors  required 
for  our  solvation  model.  ^ ’2  This  model  may  be  denoted  as  a  fragment  model,  as  it  estimates  the 
value  of  a  given  solvent  deseriptor  as  a  sum  of  eontributions  from  various  types  of  fragments  in 
the  solute.  These  fragments  are,  in  some  eases,  ambiguous,  partieularly  for  N-eontaining  solutes. 
Beeause  the  fragment  model  proposed  by  Platts  and  eoworkers  ean  be  ambiguous  and  beeause  it 
is  not  available  for  two  (the  index  of  refraetion  and  the  maeroeopie  surfaee  tension)  of  the  four 
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required  deseriptors  in  our  solvation  model,  we  are  in  the  proeess  of  developing  several  models 
eomparable  to  this  model.  In  partieular  one  model  eorrelates  solvent  deseriptor  values  to  the 
different  types  of  eovalent  bonds  in  the  solvent,  and  another  model  eorrelates  these  values  to  the 
exposed  surfaee  areas  of  the  atoms  in  the  solvent.  The  former  model  is  also  a  fragment  model, 
but  the  fragments  are  unambiguous.  In  the  exposed  surfaee  area  model,  the  only  required  input 
is  the  three  dimensional  geometry  of  the  solute.  Preliminary  results  show  that  these  models 
perform  as  well  as  or  better  than  the  fragment  model  of  Platts  and  eoworkers. 
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•  "Fast  and  Accurate  Methods  for  Predicting  Solvent  Descriptors,"  Thompson,  J.  D.;  Cramer,  C. 
J.;  Truhlar,  D.  G. 
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Method  Development  and  Simulations  of  Physical  Processes  and  Chemical  Reactions  in  the 
Various  Phases 

(Donald  L.  Thompson,  Department  of  Chemistry,  Oklahoma  State  University,  Stillwater, 
Oklahoma) 

The  foeus  of  the  researeh  at  Oklahoma  State  University  is  the  development  of  praetieal 
methods  for  simulations  and  rate  ealculations  of  physieal  and  ehemical  proeesses  in  the  gas, 
liquid,  and  solid  phases  of  energetic  materials.  The  work  for  this  report  period  focused  on  the 
following  projects: 

•  Development  of  accurate  methods  for  simulating  melting  and  predicting  melting  points  of 
organic  solids. 

•  Development  practical  methods  for  computing  reaction  rates  of  reactions  in  condensed  phases. 

•  Development  of  methods  for  fitting  ab  initio  potential  energy  surfaces  and  for  performing 
direct  dynamics  simulations. 

•  Using  quantum  chemistry  methods  to  determine  the  reaction  pathways  for  unimolecular  and 
bimolecular  reactions  of  energetic  molecules. 

Predictive  models  for  condensed  phase  processes 

A  significant  result  of  our  work  during  this  report  period  is  the  development  and 
demonstration  of  methods  for  simulating  solid-to-liquid  phase  transitions  in  energetic  materials 
to  predict  their  melting  points.  Prior  to  this  grant,  working  with  Drs.  Betsy  Rice  and  Dan 
Sorescu,  we  developed  crystal  models  for  organic  materials  that  accurately  describe  a  wide  range 
of  energetic  materials.  In  the  current  grant  we  have  used  these  models  in  studies  of  melting 
transitions  of  energetic  materials.  The  simulation  of  melting  and  the  prediction  of  melting  points 
have  been  considered  difficult  problems.  The  difficulty  stems  from  the  free  energy  barrier  for 
the  formation  of  a  liquid-solid  interface;  which  can  cause  superheating  in  a  perfect  crystal  and 
thus  an  overestimation  of  the  melting  point  of  the  real  material.  We  have  investigated  various 
ways  of  performing  the  simulations  to  determine  the  most  practical  one  for  use  in  studies  of 
energetic  materials.  Prior  to  our  work,  the  simulations  had  been  restricted  to  atomic  solids  such 
as  rare  gases  and  metals.^  We  first  carried  out  a  series  of  calculations  for  argon^  to  test  methods 
which  we  then  extended  to  nitromethane,^  RDX,"^  and  PETN.^ 

The  introduction  of  vacancies  in  a  crystal  eliminates  the  free  energy  barrier  that  causes 
superheating  in  a  perfect  crystal.  Alternatively,  the  simulations  can  be  initiated  with  the  solid  in 
contact  with  liquid.  Molecular  dynamics  studies  of  melting  of  nitromethane  using  the 
intermolecular  interaction  potential  of  Sorescu  et  al.^  have  been  carried  out  by  two  methods:  (1) 
void-nucleated  melting  with  the  gradual  heating  of  the  lattice  and  (2)  equilibration  of  coexisting 
liquid  and  solid  phases.  The  results  are  in  near  agreement  with  each  other;  the  small  difference  is 
attributed  to  the  hysteresis  effect  associated  with  the  direct  heating  process.  The  computed 
values  of  the  melting  temperature  are  in  good  agreement  with  the  experimental  data  at  various 
values  of  pressure  ranging  from  1  atm  to  30  kbar.  Since  the  void-induced  melting  approach  is 
much  easier  to  use,  we  have  adopted  it  as  the  practical  approach  for  predicting  melting  points  of 
energetic  materials.  Following  this  work,  which  used  our  well-tested  force  field  for 
nitromethane,  we  have  extended  the  studies  to  RDX  and  PETN,  which  we  believe  will  provide 
critical  tests  of  the  models  and  the  melting  approaches.  We  find  that  the  results  are  sensitive  to 
details  of  the  potential  and  that  for  large,  complex  molecules  the  simulations  are  somewhat  more 
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demanding;  however,  the  results,  which  are  currently  being  prepared  for  publication,  show  that 
approaches  that  we  are  using  will  be  effective  and  practical  for  predicting  melting  points. 

This  work  is  being  done  in  collaboration  with  Dr.  Betsy  Rice  at  ARL.  She  will 
incorporate  these  methods  into  the  general  simulation  codes  she  is  developing  for  Army 
applications. 

Practical  Rate  Calculation  Methods 

While  the  main  thrust  of  the  work  has  been  explicit  MD  simulations  of  processes  in 
condensed  phase  energetic  materials,  we  are  also  interested  in  developing  accurate  methods  for 
rate  calculations  that  complement  the  results  of  simulations,  especially  methods  that  can  be  used 
to  either  save  computer  time  or  to  treat  reactions  that  cannot  be  dealt  with  in  a  simulation.  Thus, 
we  have  explored  methods  for  computing  rates  in  the  liquid  phase.  We  have  used  relatively 
simply  reactions  and  performed  the  studies  for  reactions  in  simple  solvents;  however,  the 
approach  can  be  extended  to  more  complex  systems  of  practical  interests  in  energetic  materials 
research. 

The  solvent  effects  on  rates  of  chemical  reactions  have  been  extensively  studied  over  the 
years.  Much  of  the  theoretical  work  is  based  on  stochastic  methods  such  as  the  Langevin  theory. 
These  methods  allow  for  computations  of  rate  constants  without  explicit  consideration  of  the 
detailed  dynamics  of  the  solvent,  and  are  thus  powerful  for  treating  complex  systems  where  full¬ 
dimensional  dynamical  calculations  are  too  costly.  For  instance,  in  the  Langevin  theory,  the 
solvent  influence  on  the  reacting  particle  is  modeled  by  a  random  fluctuating  force  and  a  friction 
kernel.  Although  these  stochastic  methods  are  useful  for  qualitatively  understanding  reactions  in 
condensed  phases,  rate  constants  cannot  be  directly  obtained  for  a  given  realistic  system  because 
quantities  involved  in  the  models  concerning  the  solvent  effects  are  generally  unknown  and  must 
be  determined  by  examining  the  detailed  dynamics  of  the  system.  Molecular  dynamics  (MD) 
simulations  provide  a  useful  means  for  realistically  modeling  detailed  dynamics  and  elucidating 
reaction  mechanisms.  However,  a  full-dimensional  MD  study  is  computationally  expensive  and 
may  not  be  feasible  for  complex  systems.  An  alternative  approach  is  to  map  the  real  system  to  a 
stochastic  model  by  performing  a  few  MD  calculations  to  determine  the  fitting  parameters  in  the 
model.  We  have  explored  such  an  approach. 

The  basic  idea  is  that  the  reacting  molecule  is  subjected  to  random  collisions  with  a  mean 
frequency  a.  Between  collisions,  the  dynamics  of  the  molecule  is  governed  solely  by  the 
equations  of  motion  for  the  molecule  itself  The  effects  of  the  solvent  are  modeled  by  random 
collisions.  Under  the  strong  collision  assumption,  the  velocities  of  the  molecule  are  randomized 
after  each  collision  by  sampling  from  the  Boltzmann  distribution  while  the  coordinates  are  held 
fixed.  Since  only  the  motion  of  the  molecule  is  explicitly  treated,  the  method  is  far  less  costly 
than  full-dimensional  MD  simulations.  The  parameter  used  to  model  the  solvent  effects  is  the 
collision  frequency  a,  which  is  related  to  the  liquid  density  p  of  the  system.  There  have  been 
studies  where  both  MD  simulations  and  the  stochastic  dynamics  method  were  employed  on  the 
same  system;  but  no  attempt,  to  our  knowledge,  had  previously  been  made  to  quantitatively 
relate  a  and  p  for  any  system. 

We  first  investigated  this  idea  with  an  application  to  cis-trans  isomerization  in  HONO, 
which  is  one  of  the  simplest  molecules  that  exhibit  cis-trans  isomerization  and  one  that  has  been 
extensively  studied  experimentally  and  theoretically.  We  used  a  model  system  with  a  HONO 
molecule  embedded  in  liquid  krypton  and  performed  MD  simulations  to  obtain  isomerization 
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rates  at  several  liquid  densities.  We  also  eomputed  rate  eonstants  using  the  stoehastic  dynamies 
method  for  a  wide  range  of  eollision  frequeneies  a  for  these  two  sets  of  eomputed  rates.  We 
then  examine  the  relationship  between  the  liquid  density  p  and  the  collision  frequency  a. 
Comparisons  of  the  two  sets  of  the  computed  rates  show  that  for  a  wide  range  of  liquid  densities 
there  is  a  simple  linear  relation  between  the  liquid  density  p  and  the  collision  frequency  a,  that 
is,  a  =  cp.  This  suggests  that  once  the  constant  c  is  determined  from  a  molecular  dynamics 
calculation  at  a  single  density,  the  reaction  rates  can  be  obtained  from  stochastic  dynamics 
calculations  for  the  entire  range  of  liquid  densities  where  a  =  c/?  holds.  The  applicability  of  the 
combined  molecular  dynamics  and  stochastic  dynamics  approach  provides  a  practical  means  for 
obtaining  rate  constants  at  considerable  savings  of  computer  time  compared  to  that  required  by 
using  full-dimensional  molecular  dynamics  simulations  alone.  The  results  of  this  work  have 
been  published.^  Given  the  success  of  this  study,  we  next  carried  out  a  similar  study  of  bond 
fission  in  MONO  (i.e.,  MONO  ^  OH  +  NO). 

Intuitively,  the  solvent  effects  on  dissociation  reactions  should  be  different  from  those  for 
isomerizations;  the  cage  effect  due  to  the  solvent  that  forces  recombination  should  be  a 
prominent  factor  in  dissociation  reactions.  Thus,  the  relation  between  a  and  p  as  well  as  the 
dependence  of  the  reaction  rate  on  p  are  expected  to  be  different  from  that  for  isomerization.  To 

establish  the  relation  between  a  and  p,  we  have  performed  two  sets  of  calculations.  We 
calculated  the  rate  as  a  function  of  a  using  the  stochastic  dynamics  method.  We  also  calculated 
the  rate  at  four  liquid  densities  by  using  the  full-dimensional  MD  simulations  with  a  HONO 
molecule  inside  liquid  Kr.  The  relation  between  a  and  p  was  found  by  pairing  a  and  p  that 
have  the  same  rate.  In  our  study  of  the  cis-trans  isomerization  of  HONO  in  liquid  Kr,  we  found 
that  a  and  p  obey  a  simple  linear  relation,  a  =  cp  .  For  the  0-N  bond  dissociation  of  HONO  in 
liquid  Kr,  we  found  that  the  relationship  is  more  complicated,  but  is  accurately  described  by  an 
analytical  expression.  We  have  also  found  that  the  solvent  effects  are  different  for  the 
isomerization  and  dissociation  reactions.  For  isomerization  the  rate  is  in  the  low-collision  regime 
and  increases  with  increasing  liquid  density,  whereas  for  dissociation  the  rate  is  in  the  high- 
collision  regime  and  decreases  with  increasing  liquid  density.  A  manuscript  has  been  submitted 
for  publication. 

In  both  cases,  based  on  one  MD  simulation  at  a  single  density,  the  reaction  rates  can  be 
obtained  from  stochastic  dynamics  for  a  wide  range  of  liquid  densities.  The  approach  thus 
provides  tremendous  savings  of  computational  time  compared  to  full-dimensional  MD 
simulations.  The  results  to  date  show  that  the  method  is  applicable  to  isomerization  and 
dissociation  reactions  of  small  molecules  in  rare-gas  liquids.  In  future  work  we  will  extend  the 
studies  to  larger,  more  complicated  systems,  including  more  realistic  energetic  molecules  and  in 
reactions  in  neat  liquids. 

Reaction  Pathways  for  Unimolecular  Dissociations  of  Energetic  Molecules 

A  major  challenge  to  theoretical  chemists  and  a  particular  need  in  energetic  materials 
research  is  the  development  of  methods  for  simulations  and  rate  calculations  for  sequential, 
branching  reactions.  The  practical  problems  range  from  soot  formation  to  the  decomposition  of 
an  energetic  material  such  as  RDX.  We  have  long  been  interested  in  this  problem  within  the 
context  of  unimolecular  decomposition  reactions.  We  are  pursuing  two  lines  of  research  that  we 
expect  to  join  when  both  are  more  mature.  First,  we  are  using  quantum  chemistry  to  explore  the 
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transition  states  and  intermediates  in  the  ehemieal  deeomposition  of  energetie  moleeules.  The 
other  aspeet  of  this  work  is  purely  developmental  at  this  stage  and  foeuses  on  methods  for  fitting 
ab  initio  potential  energy  surfaees  and  for  performing  direet  dynamies  simulations  (a  brief 
mention  of  this  work  is  given  below). 

During  this  report  period  quantum  ehemistry  studies  of  the  deeomposition  pathways  of 
1,3,3-trinitroazetidine  (TNAZ)  and  dimethylnitramine  have  been  carried  out.  Possible 
decomposition  transition-states,  intermediates,  and  products  of  TNAZ  were  identified  and  the 
structures,  energies,  and  vibrational  frequencies  were  determined  at  the  B3LYP/6-31G(<i,j!7)  level. 
Four  major  pathways  are  apparent.  Two  pathways  are  initiated  by  the  fission  of  the  N-NO2  and 
C-NO2  bonds  to  yield  radical  intermediates,  while  the  other  two  pathways  involve  the  molecular 
elimination  of  HONO.  Energy  profiles  for  the  pathways  and  possible  routes  to  some  of  the 
experimentally  observed  species  of  TNAZ  decomposition  have  been  determined.  The  energies 
required  to  initiate  the  NO2  bond  fission  pathways  are  4  to  8  kcal/mol  lower  than  the  HONO 
elimination  pathways.  In  the  gas  phase,  the  NO2  elimination  pathways  will  be  the  dominant 
routes  for  TNAZ  decomposition.  In  the  condensed  phase,  however,  this  trend  may  be  reversed. 
This  work  has  been  published:  S.  Alavi,  L.  M.  Reilly,  and  D.  L.  Thompson,  J.  Chem.  Phys.  119, 
8297  (2003). 

Examination  of  the  reported  experimental  and  theoretical  results  from  several  researchers 
shows  that  there  is  an  significant  dispersion  in  the  values  of  the  Arrhenius  parameters  of  the  rate 
constants  for  the  unimolecular  decomposition  of  DMNA.  The  structures  and  energies  of  the 
reactant,  intermediates,  products,  and  transition  states  of  the  initial  stages  of  DMNA 
decomposition  of  have  been  determined  by  quantum  chemical  calculations  at  several  levels  of 
theory.  Kinetic  parameters  for  these  decomposition  steps  have  been  obtained  using  the  RRKM 
formalism  for  the  range  300-5000  K.  The  pathways  considered  are  NO2  elimination,  HONO 
elimination,  and  NNO2-NONO  rearrangement.  The  NO2  elimination  channel  is  the  main  channel 
of  gas  phase  decomposition  of  DMNA  in  the  300-5000  K  range  of  temperatures.  The  HONO 
elimination  channel  has  the  next  lowest  decomposition  energy  but  at  high  temperatures  the  nitro- 
nitrite  rearrangement  competes  with  it. 

The  rate  constants  as  functions  of  temperature  for  NO2  elimination,  HONO  elimination, 
and  nitro-nitrite  arrangement  have  been  obtained  from  values  of  the  energies  and  frequencies 
using  RRKM  theory;  Arrhenius  plots  are  shown  in  the  figure  below,  where  the  curves  are 
identified  as  follows:  (— )  NO2  elimination,  ( — )  HONO  elimination,  and  (•••)  nitro-nitrite 
rearrangement. 
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The  Arrhenius  parameters  for  NO2  elimination  are  in  good  agreement  with  the  values  of  previous 
theoretieal  studies.  The  HONO  elimination  ehannel  has  a  lower  value  of  the  pre-exponential 
faetor  eompared  to  NO2  elimination  beeause  the  transition  state  for  the  HONO  elimination 
pathway  is  entropieally  unfavorable.  Clearly,  the  NO2  elimination  is  the  main  deeomposition 
ehannel  of  DMNA  in  the  range  of  300-5000  K.  The  HONO  elimination  has  a  higher  rate  than 
the  nitro-to-nitrite  rearrangement  but  at  high  temperatures  the  rates  of  these  two  reaetions 
beeomes  eomparable.  The  experimental  evidenee  is  that  the  presenee  of  NO  in  the  DMNA 
deeomposition  produets  is  due  to  the  nitro-nitrite  rearrangement  followed  by  breaking  of  the 
weak  NO-NO  bond.  Our  results  show  that  at  low  temperatures,  this  ehannel  is  the  third  in 
importanee  and  beeomes  eomparable  with  the  HONO  elimination  ehannel  only  at  higher 
temperatures. 

An  important  part  of  predieting  the  behavior  of  energetie  materials  is  better  praetieal 
methods  for  eomputing  the  rates  of  ehemieal  reaetions.  A  eritieal  problem  is  the  formulation  of 
PESs.  For  purely  predietive  methods  it  neeessary  that  the  PESs  be  based  on  quantum  ehemistry 
ealeulations,  but  this  remains  a  major  obstaele  to  applieations  to  large  polyatomie  moleeules  and 
radieals.  The  most  straightforward  way  to  make  use  of  ab  initio  results  is  to  simply  eompute  the 
forees  “on  the  fly”  during  a  simulation;  however,  the  diffieulty  here  is  that  studies  are  restrieted 
to  the  lower-level  quantum  methods,  whieh  are  often  not  suffieiently  aeeurate  for  reliable 
predietions  of  rates.  The  eommon  alternative  to  “direet  dynamies”  is  to  fit  the  ab  initio  energies 
to  an  analytieal  funetion,  whieh  allows  for  mueh  more  rapid  evaluations  of  the  forees  during 
simulations  as  well  as  sealing  of  the  ab  initio  results  to  eorreet  for  the  inaeeuraeies  inherent  in  a 
low-level  theory.  Fitting  a  global  PES  is  extremely  tedious  and  not  readily  generalizable,  thus  a 
huge  investment  of  human  labor  is  required  in  eaeh  ease.  Furthermore,  the  analytieal  forms  are 
not  suffieiently  flexible  to  aeeurately  fit  the  ab  initio  points.  This  was  our  motivation  for 
suggesting,  in  the  1970’s,  that  loeal  fitting  funetions  should  be  used  to  fit  ab  initio  energies.  We 
proposed  using  eubie  splines,  whieh  provide  the  desired  flexibility  to  yield  aeeurate,  smooth  first 
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and  seeond  derivatives  but  requires  a  fairly  dense  grid  of  ab  initio  points.  The  modified 
Shepard  methods  suggested  by  Isehtwan  and  Collins^  in  1994  are  superior  to  eubie  splines.  They 
have  eombined  surfaee  fitting  with  MD  simulations  in  an  iterative  seheme  for  sueeessively  and 
automatieally  improving  the  PES.  This  proeedure  seleets  the  loeations  for  additional  ab  initio 
ealeulations  in  regions  of  eonfiguration  spaee  that  are  dynamieally  important.  The  proeedure  is 
simple  and  readily  automated  but  gradients  and  Hessians,  whieh  are  not  readily  available  in  the 
high-level  ab  initio,  are  required  at  every  point.  We  have  developed  interpolative  moving  least- 
squares  (IMLS)  fitting  proeedures  that  promise  to  be  more  useful  than  other  methods,  mainly 
beeause  it  does  not  require  derivatives. 

During  this  report  period  we  began  applieations  of  the  IMLS  methods  to  reaetions 
important  in  the  deeomposition  of  nitramines.  We  have  begun  with  a  study  of  the  unimoleeular 
deeomposition  of  H2CN,  whieh  is  eurrently  under  study  in  Professor  Paul  Dagdigian’s  laboratory 
at  Johns  Hopkins.  We  are  doing  this  work  in  eollaboration  with  Drs.  Larry  Harding  and  A1 
Wagner  at  ANL.  Ab  initio  energies  have  been  fit  and  rates  of  dissoeiation  are  being  earried  out. 
These  results  are  being  generated  in  advanee  of  Professor  Dagdigian’s  measurements  of  the 
rates,  thus  they  will  be  truly  predietive.  The  immediate  goal  is  to  establish  the  aeeuraey  of  this 
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approach  for  a  significant  chemical  reaction.  We  plan  to  extend  the  work  to  include  bimolecular 
reaetions  of  H2CN  with  other  radieals  and  small  moleeules. 

Our  goal  is  to  develop  software  that  manages  the  quantum  ehemistry  caleulations  in 
direet  response  to  the  needs  of  a  MD  simulation  to  provide  an  ab  initio  prediction  of  a  specified 
reaetion  rate.  If  sueeessful,  and  the  prospeets  are  good,  this  eould  make  MD  studies  of  reaetions 
a  praetieal  tool  for  all  who  are  interested  in  reaetion  kineties  rather  than  a  teehnique  limited  to 
use  by  experts. 
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Plans:  The  following  will  be  the  focus  of  our  group  for  the  coming  year: 

•  Continuation  of  studies  of  melting  of  energetie  materials 

•  Further  development  of  methods  for  rate  ealeulations  for  eondensed-phase  reaetions 

•  Refinement  and  applieations  of  methods  for  fitting  ab  initio  PESs  for  reaetions 

•  Rate  calculations  and  dynamics  calculations  for  the  chemical  decomposition  of  energetic 
molecules 

•  Simulations  of  shoeked  solids  and  liquids 

Personnel: 

Ms.  Lisa  M.  Riley  (Graduate  Student,  Supported  by  MURI  funds) 

Dr.  Paras  M.  Agrawal  (Researeh  Assoeiate,  Supported  by  MURI  funds.) 

Dr.  Saman  Alavi  (Research  Associate,  Partially  supported  by  MURI  funds.) 

Dr.  Yin  Guo  (Research  Associate,  Partially  supported  by  MURI  funds.) 
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the  Melting  of  Nitromethanef  J.  Chem.  Phys.  119,  9617  (2003). 

•  Saman  Alavi,  Lisa  M.  Reilly,  and  Donald  L.  Thompson,  ""Theoretical  Predictions  of  the 
Decomposition  Pathways  of  1 ,3 ,3-Trinitroazetidine  (TNAZ),"  J.  Chem.  Phys.  119,  8297 
(2003). 

•  Yin  Guo  and  Donald  L.  Thompson,  ""On  Combining  Molecular  Dynamics  and  Stochastic 
Dynamics  Simulations  to  Compute  Reaction  Rates  in  Liquids,"  J.  Chem.  Phys.  120,  898 

(2004). 
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2004,  in  press. 
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Non  Peer  Reviewed  Papers:  (4) 

Book  Chapters: 

•  '^Structure  and  Density  Predictions  for  Energetic  Materials  f  J.  R.  Holden,  Z.  Du  and  H.  L. 
Ammon,  in  Energetic  Materials,  Part  1:  Decomposition,  Crystal  and  Molecular  Properties, 
P.  Politzer  and  J.  S.  Murray,  eds.,  pp.  183-213,  Elsevier,  Amsterdam,  2003. 

Papers  Published  in  Conference  Proceedings: 
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J.;  Truhlar,  D.  G. 
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•  R.  M.  Martin,  ""Electronic  Structure:  Basic  Theory  and  Practical  Methods",  Cambridge 
University  Press,  2004 

Faculty  Supported:  (9) 

John  Adams 
Herman  Ammon 
Rodney  Bartlett 
Donald  W.  Brenner 
David  Ceperley 
Chistopher  J.  Cramer 
Richard  Martin 
Donald  L.  Thompson 
Donald  G.  Truhlar 


33 


Graduate  Students  Supported:  (14) 

Ed  Wells 
William  Mattson 
Edward  Bukhman 
Elsa  M.  Riley 
DeCarlos  Taylor 
Tom  Henderson 
Anthony  Yau 
Josh  MeClellan 
Tom  Hughes 
Euis  Galiano 
Staeie  Gregory 
All  Siavosh-Haghighi 
Jason  D.  Thompson 
Tomas  Szabo 

Postdoctoral  Fellows  Supported:  (5) 

Zuyue  Du 
Simone  Chiesa 
Paras  M.  Agrawal 
Saman  Alavi 
Yin  Guo 

Other  Research  Staff:  (2) 

Sayta  Prasad,  sabbatical  associate  (Physics,  Ranchi  Elniv,  India) 
Xinlu  Chen  (visiting  professor) 

Kai  Wang  (senior  researcher) 

Undergraduate  Students  Supported:  (1) 

Nicole  Dueker 

PhDs  Awarded: 

W.  D.  Mattson,  thesis.  University  of  Illinois,  2003. 

Substantial  Interactions  with  DoD  and  Other  Researchers: 

Dr.  Betsy  Rice  (ARE) 

Dr.  Ed  Byrd  (ARE) 

Dr. William  Mattson  (ARE) 

Dr.  Alfred  Stern  (NAVSEA,  Indian  Head,  MD) 

Dr.  Wiliam  Koppes  (NAVSEA,  Indian  Head,  MD) 

Dr.  Rao  Surapaneni  (Picatinny  Arsenal) 

Dr.  Paritosh  Dave  (Picatinny  Arsenal) 

Dr.  Robert  Chapman  (NAWC,  China  Eake,  CA) 

Dr.  Philip  Eaton  (Univ.  of  Chicago) 

Dr.  Harold  Shechter  (The  Ohio  State  Univ.) 

Dr.  Steve  Bunte  (ARE) 

Dr.  Shashi  Kama  (ARE) 

Dr.  Cary  Chabalowski  (ARE) 


34 


